This notebook is the relative cumulative frequency analysis of metabolic cage data.
Data was generated previously when I ran PCA, see notebook .
The idea behind this was to apply a method to study large metabolic datasets quantitatively. Based on the Gordon paper, we can use relative cumulative frequency to find frequency of data binned by VO2 or RQ reading levels, then calculate EC50. We can use EC50 values for every mouse in a simple t-test or ANCOVA to compare groups. In other words, this can be thought of as a kind of weighted average.
There are some assumptions that need to be considered. First are the flaws in the dataset — it is difficult to draw any conclusions when the datasets are flawed to begin with. There are 12 SPF mice, 5 of which are WT and 7 of which are L-Bmal1-KO. There are 6 GF mice, of which mouse #4 was contaminated and cannot be used, and of the remaining 5, 2 are WT and 3 are L-Bmal1-KO. The mix of genotypes prevents any far-reaching conclusions. Furthermore, due to limitations in the tools I am using (there needs to be the same number of samples for both groups), I am stuck to using only 5 SPF mice, to compare with 5 GF mice. Overcoming this requirement (i.e. having unequal groups) will be a good idea in the future.
Thus, this notebook serves as a proof-of-concept, that relative cumulative frequency can be applied, EC50 calculated, and statistics run.
Relative cumulative frequency
See the paper for more information . This is a method where data is first ordered by increasing value, then checked for increases based on some increment, e.g. for every increment with an increase, a frequency is recorded. Finally, frequency is cumulated. This is all accomplished by the function rcf. In order to run this, however, the range of data needs to be checked. find_range can be run to check the range on a vector, and will return a recommended range and increment value to submit to rcf.
EC50
By fitting frequencies to intervals, a sigmoidal regression is found. This is dependent on a normal distribution — that is, the mouse should have less frequencies towards the extremes, and the most expression in the middle. For RQ and VO2, this is true, the mouse will not spend as much time on the extremes of the values. For other data types, the nature and distribution of the data should be considered before using in this analysis.
t-test and ANCOVA
For the Welch t-test, I run an F-test first (var.test), which should be not significant.
These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).
| Start of Met Cage run |
Metabolic Cage # |
Ear Tag | DOB | Gender | Genotype | Body Mass on BMR day (g) |
|---|---|---|---|---|---|---|
| 3/10/2020 | 1 | 1555 | 11/18/19 | M | WT | 26.4 |
| 3/10/2020 | 2 | 1692 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 3 | 1693 | 12/06/19 | M | WT | 24.7 |
| 3/10/2020 | 4 | 1695 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 5 | 1699 | 12/06/19 | M | WT | 28.7 |
| 3/10/2020 | 6 | 1556 | 11/18/19 | M | L-Bmal1-KO | 26.2 |
| 3/10/2020 | 9 | 1557 | 11/18/19 | M | L-Bmal1-KO | 27.1 |
| 3/10/2020 | 10 | 1558 | 11/18/19 | M | L-Bmal1-KO | 26.9 |
| 3/10/2020 | 11 | 1559 | 11/18/19 | M | L-Bmal1-KO | 24.7 |
| 3/10/2020 | 12 | 1691 | 12/06/19 | M | L-Bmal1-KO | 28 |
| 3/10/2020 | 13 | 1694 | 12/06/19 | M | L-Bmal1-KO | 26.4 |
| 3/10/2020 | 14 | 1700 | 12/06/19 | M | L-Bmal1-KO | 27.7 |
| Cage | Ear Tag | Genotype | GF condition, day 6 | GF condition, day 10 |
|---|---|---|---|---|
| 1 | 5650 | L-Bmal1-KO | GF | contaminated |
| 2 | 976 | WT | GF | GF |
| 3 | 977 | WT | GF | contaminated |
| 4 | 990 | L-Bmal1-KO | likely contaminated | likely contaminated |
| 5 | 979 | L-Bmal1-KO | GF | contaminated |
| 6 | 978 | L-Bmal1-KO | GF | contaminated |
For the first half of the experiment, only cage 4 had contamination. For the second leg of the experiment, cages 1,3,5 and 6 were contaminated, and it is likely cage 4 was also contaminated.
See the 16S rDNA gel below.
reqpkg <- c("docstring","magrittr","dplyr","ggplot2","ggpubr","reshape2","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "docstring"
## [1] '1.0.0'
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "reshape2"
## [1] '1.4.3'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
theme_set(theme_pubr())
signif.floor <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- floor(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
signif.ceiling <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- ceiling(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
find_range <- function(data, channel, segments=10) {
data <- data %>% dplyr::select(contains(channel)) %>% melt()
r <- range(data$value)
cut <- max(data$value)-min(data$value)
out <- list("range"=r, "cut"=cut/segments)
return(out)
}
rcf <- function(data, channel, id, min, max, b=0.1) {
#' Relative cumulative frequency function
#'
#' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
#'
#' @param data numeric vector
#' @param channel character. Channel to label the rcf output.
#' @param b double. Number to cut vector.
#' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
breaks <- seq(min, max, by=b)
# breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
duration.cut <- cut(data, breaks, right = F)
duration.freq <- table(duration.cut)
duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
duration.cumrelfreq = duration.cumfreq/length(data)
out <- duration.cumrelfreq %>% as.data.frame()
out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
out$id <- rep(id, nrow(out))
return(out)
}
EC50 <- function(x) {
out <- lapply(x, function(y) {
nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
}) %>%
sapply(function(z) {
getEstimates(z, 0.5)$x
})
out %<>% set_rownames(NULL)
return(out)
}
lm_eqn <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
lm_eqn_str <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
return(eq)
}
remove_outliers <- function(x, na.rm = TRUE, ...) {
qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
H <- 1.5 * IQR(x, na.rm = na.rm)
y <- x
y[x < (qnt[1] - H)] <- NA
y[x > (qnt[2] + H)] <- NA
y
}
ANCOVA <- function(df) {
# Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
cat("linear fit of dependent variable with covariates\n")
p1 <- ggplot(df, aes(df[[2]], df[[3]])) +
geom_smooth(method = "lm", se = FALSE) +
geom_point() +
# geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
ggtitle("Linear model - DV and Covar") +
xlab("DV") + ylab("Covar")
print(p1)
m <- lm(df[[3]] ~ df[[2]], df)
summary(m) %>% print()
lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
plot1 <- ggplot(df, aes(x=df[[1]], y=df[[2]])) +
geom_boxplot() +
# ggtitle("DV vs IV") +
xlab("IV") + ylab("DV")
plot2 <- ggplot(df, aes(x=df[[1]], y=df[[3]])) +
geom_boxplot() +
# ggtitle("Covar vs IV") +
xlab("IV") + ylab("Covar")
p2 <- ggarrange(plot1, plot2,
labels = c("DV vs IV", "Covar vs IV"),
ncol = 2)
# grid.arrange(plot1, plot2, ncol=2)
print(p2)
cat("\n_____________________\n")
cat("Some stats on the variables\n")
cat("DV - IV:\n")
by(df[[2]], df[[1]], stat.desc) %>% print()
cat("Covar - IV:\n")
by(df[[3]], df[[1]], stat.desc) %>% print()
cat("levene's test\n")
leveneTest(df[[2]], df[[1]], center = median) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV\n")
model <- aov(df[[2]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of covariate with IV, independent of DV\n")
model <- aov(df[[3]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV, normalized to covariate\n")
c <- df[[2]]/df[[3]]
model <- aov(c ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANCOVA\n")
model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
summary(model) %>% print()
out <- Anova(model, type= "III")
print(out)
return(out)
}
check_ANCOVA <- function(res) {
if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
} else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
} else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
}
}
SPF_dark <- read.csv2("data/R/SPF_dark.csv", sep = ",") %>%
set_colnames(paste0("SPF_",colnames(.)))
SPF_light <- read.csv2("data/R/SPF_light.csv", sep = ",") %>%
set_colnames(paste0("SPF_",colnames(.)))
SPF_covar <- read.csv2("data/R/SPF_covar.csv", sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
GF_dark <- read.csv2("data/R/GF_dark.csv", sep = ",") %>%
dplyr::select(-contains(c("_M_4","_R_4","_Mnz_4"))) %>%
set_colnames(paste0("GF_",colnames(.)))
GF_light <- read.csv2("data/R/GF_light.csv", sep = ",") %>%
dplyr::select(-contains(c("_M_4","_R_4","_Mnz_4"))) %>%
set_colnames(paste0("GF_",colnames(.)))
GF_covar <- read.csv2("data/R/GF_covar.csv", sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
light <- SPF_light %>%
bind_cols(GF_light) %>%
dplyr::select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark <- SPF_dark %>%
bind_cols(GF_dark) %>%
dplyr::select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total <- bind_rows(light, dark)
Select a tab above to view
total.VO2 <- total %>% dplyr::select(contains("VO2"))
total.VO2.SPF <- total.VO2 %>% dplyr::select(contains("SPF"))
total.VO2.GF <- total.VO2 %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(total.VO2.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
total.VO2.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(total.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(total.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(tmp){
nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in complete diurnal cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.68268673163679"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.54667034693964"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.4691, num df = 4, denom df = 4, p-value = 0.7184
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1529618 14.1102621
## sample estimates:
## ratio of variances
## 1.469126
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = 2.1814, df = 7.7213, p-value = 0.06193
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.008805302 0.284815805
## sample estimates:
## mean of x mean of y
## 1.680105 1.542100
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7915 -0.5634 -0.0086 0.3746 3.6986
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.7421 7.8325 3.414 0.00917 **
## df[[2]] 0.4766 4.8496 0.098 0.92413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.733 on 8 degrees of freedom
## Multiple R-squared: 0.001206, Adjusted R-squared: -0.1236
## F-statistic: 0.009659 on 1 and 8 DF, p-value: 0.9241
##
## y = 27 + 0.48x, r^2 = 0.00121
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.000000000 0.000000000 0.000000000 1.394711940 1.611752055 0.217040115
## sum median mean SE.mean CI.mean.0.95 var
## 7.710499431 1.593015270 1.542099886 0.040261519 0.111783896 0.008104949
## std.dev coef.var
## 0.090027493 0.058379806
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 1.57230931 1.86329413 0.29098482
## sum median mean SE.mean CI.mean.0.95 var
## 8.40052569 1.65195795 1.68010514 0.04879998 0.13549046 0.01190719
## std.dev coef.var
## 0.10912007 0.06494836
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0129 0.9123
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.04761 0.04761 4.758 0.0607 .
## Residuals 8 0.08005 0.01001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 1.698e-04 1.698e-04 15.31 0.00446 **
## Residuals 8 8.868e-05 1.109e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.04761 0.04761 5.491 0.0516 .
## df[[3]] 1 0.01935 0.01935 2.232 0.1788
## Residuals 7 0.06070 0.00867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.008611 1 0.9931 0.35217
## df[[1]] 0.066812 1 7.7053 0.02746 *
## df[[3]] 0.019352 1 2.2319 0.17883
## Residuals 0.060696 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0274620158378014"
light.VO2 <- light %>% dplyr::select(contains("VO2"))
light.VO2.SPF <- light.VO2 %>% dplyr::select(contains("SPF"))
light.VO2.GF <- light.VO2 %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(light.VO2.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
light.VO2.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(light.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(light.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(tmp){
nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.4791478793743"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.35106555611578"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.6667, num df = 4, denom df = 4, p-value = 0.6328
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1735299 16.0076157
## sample estimates:
## ratio of variances
## 1.666674
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = 2.1057, df = 7.5294, p-value = 0.07048
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.01294961 0.25464656
## sample estimates:
## mean of x mean of y
## 1.472360 1.351512
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8906 -0.7058 -0.2072 0.6747 3.4101
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.480 7.540 4.175 0.0031 **
## df[[2]] -2.812 5.327 -0.528 0.6119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.704 on 8 degrees of freedom
## Multiple R-squared: 0.03366, Adjusted R-squared: -0.08713
## F-statistic: 0.2787 on 1 and 8 DF, p-value: 0.6119
##
## y = 31 + -2.8x, r^2 = 0.0337
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.000000000 0.000000000 0.000000000 1.241538599 1.450636943 0.209098344
## sum median mean SE.mean CI.mean.0.95 var
## 6.757558245 1.375201548 1.351511649 0.035144752 0.097577474 0.006175768
## std.dev coef.var
## 0.078586054 0.058146783
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 1.36405635 1.61172432 0.24766797
## sum median mean SE.mean CI.mean.0.95 var
## 7.36180063 1.48720545 1.47236013 0.04537177 0.12597224 0.01029299
## std.dev coef.var
## 0.10145437 0.06890595
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.344 0.5737
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.03651 0.03651 4.434 0.0683 .
## Residuals 8 0.06588 0.00823
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0001274 1.274e-04 9.084 0.0167 *
## Residuals 8 0.0001122 1.402e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.03651 0.03651 3.985 0.0861 .
## df[[3]] 1 0.00174 0.00174 0.189 0.6765
## Residuals 7 0.06414 0.00916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.026050 1 2.8431 0.13563
## df[[1]] 0.034800 1 3.7980 0.09232 .
## df[[3]] 0.001736 1 0.1895 0.67647
## Residuals 0.064139 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
dark.VO2 <- dark %>% dplyr::select(contains("VO2"))
dark.VO2.SPF <- dark.VO2 %>% dplyr::select(contains("SPF"))
dark.VO2.GF <- dark.VO2 %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(dark.VO2.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
dark.VO2.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(dark.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(dark.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(tmp){
nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
abline(h=0.5)
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.90041457533305"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.73951498204456"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.1638, num df = 4, denom df = 4, p-value = 0.8867
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1211708 11.1776467
## sample estimates:
## ratio of variances
## 1.163789
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = 2.1973, df = 7.9544, p-value = 0.05943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.009573524 0.388474314
## sample estimates:
## mean of x mean of y
## 1.919810 1.730359
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7701 -0.5941 -0.0510 0.4014 3.6662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.8554 6.4773 3.992 0.004 **
## df[[2]] 0.9066 3.5364 0.256 0.804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.727 on 8 degrees of freedom
## Multiple R-squared: 0.008148, Adjusted R-squared: -0.1158
## F-statistic: 0.06572 on 1 and 8 DF, p-value: 0.8041
##
## y = 26 + 0.91x, r^2 = 0.00815
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 1.55076022 1.85131885 0.30055863
## sum median mean SE.mean CI.mean.0.95 var
## 8.65179656 1.80215775 1.73035931 0.05861441 0.16273970 0.01717825
## std.dev coef.var
## 0.13106581 0.07574485
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 1.78102178 2.14647377 0.36545199
## sum median mean SE.mean CI.mean.0.95 var
## 9.59904853 1.90591446 1.91980971 0.06323268 0.17556205 0.01999186
## std.dev coef.var
## 0.14139256 0.07364926
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2e-04 0.9902
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.08973 0.08973 4.828 0.0592 .
## Residuals 8 0.14868 0.01859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0002812 2.812e-04 17.06 0.0033 **
## Residuals 8 0.0001318 1.648e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.08973 0.08973 6.323 0.0401 *
## df[[3]] 1 0.04934 0.04934 3.477 0.1045
## Residuals 7 0.09934 0.01419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.001463 1 0.1031 0.75754
## df[[1]] 0.137125 1 9.6624 0.01712 *
## df[[3]] 0.049339 1 3.4766 0.10451
## Residuals 0.099341 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0171186600523453"
Select a tab above to view
total.RQ <- total %>% dplyr::select(contains("RQ"))
total.RQ.SPF <- total.RQ %>% dplyr::select(contains("SPF"))
total.RQ.GF <- total.RQ %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(total.RQ.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
total.RQ.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(total.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(total.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in complete diurnal cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.840277055484309"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.846032542612996"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 0.27142, num df = 4, denom df = 4, p-value = 0.2345
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.02826007 2.60690625
## sample estimates:
## ratio of variances
## 0.2714247
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = 0.047799, df = 6.0224, p-value = 0.9634
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02793109 0.02904508
## sample estimates:
## mean of x mean of y
## 0.8444767 0.8439197
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.53056 -0.25067 -0.04439 0.39201 1.76874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.49 20.60 3.955 0.0042 **
## df[[2]] -63.95 24.40 -2.621 0.0306 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.272 on 8 degrees of freedom
## Multiple R-squared: 0.4619, Adjusted R-squared: 0.3947
## F-statistic: 6.868 on 1 and 8 DF, p-value: 0.03062
##
## y = 81 + -64x, r^2 = 0.462
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.0000000000 0.0000000000 0.0000000000 0.8098007597 0.8735776070 0.0637768474
## sum median mean SE.mean CI.mean.0.95 var
## 4.2195986627 0.8471038023 0.8439197325 0.0103345386 0.0286932792 0.0005340134
## std.dev coef.var
## 0.0231087309 0.0273826171
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.0000000000 0.0000000000 0.0000000000 0.8274327561 0.8574982941 0.0300655380
## sum median mean SE.mean CI.mean.0.95 var
## 4.2223836496 0.8485680324 0.8444767299 0.0053841326 0.0149487485 0.0001449444
## std.dev coef.var
## 0.0120392864 0.0142565046
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.6428 0.4459
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0000008 0.0000008 0.002 0.963
## Residuals 8 0.0027158 0.0003395
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 7.240e-06 7.236e-06 1.579 0.244
## Residuals 8 3.667e-05 4.584e-06
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0000008 0.0000008 0.005 0.9457
## df[[3]] 1 0.0016246 0.0016246 10.422 0.0145 *
## Residuals 7 0.0010912 0.0001559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.0278748 1 178.8168 3.065e-06 ***
## df[[1]] 0.0003705 1 2.3769 0.16704
## df[[3]] 0.0016246 1 10.4221 0.01449 *
## Residuals 0.0010912 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.0144853254058274"
light.RQ <- light %>% dplyr::select(contains("RQ"))
light.RQ.SPF <- light.RQ %>% dplyr::select(contains("SPF"))
light.RQ.GF <- light.RQ %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(light.RQ.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
light.RQ.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(light.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(light.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.786092284744787"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.789466997689693"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 0.19478, num df = 4, denom df = 4, p-value = 0.1421
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.020280 1.870769
## sample estimates:
## ratio of variances
## 0.1947798
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -0.10582, df = 5.5013, p-value = 0.9195
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.04945103 0.04543731
## sample estimates:
## mean of x mean of y
## 0.7849262 0.7869331
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.89818 -0.60505 -0.06311 1.03107 1.81545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.56 13.40 3.923 0.0044 **
## df[[2]] -31.88 17.04 -1.871 0.0983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.446 on 8 degrees of freedom
## Multiple R-squared: 0.3044, Adjusted R-squared: 0.2174
## F-statistic: 3.501 on 1 and 8 DF, p-value: 0.09825
##
## y = 53 + -32x, r^2 = 0.304
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.000000000 0.000000000 0.000000000 0.727127187 0.833409714 0.106282527
## sum median mean SE.mean CI.mean.0.95 var
## 3.934665528 0.791403096 0.786933106 0.017350154 0.048171749 0.001505139
## std.dev coef.var
## 0.038796123 0.049300408
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.0000000000 0.0000000000 0.0000000000 0.7571566743 0.8014503624 0.0442936881
## sum median mean SE.mean CI.mean.0.95 var
## 3.9246312361 0.7869148767 0.7849262472 0.0076572940 0.0212600564 0.0002931708
## std.dev coef.var
## 0.0171222298 0.0218138072
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.1621 0.3125
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000010 0.0000101 0.011 0.918
## Residuals 8 0.007193 0.0008992
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.390e-06 5.393e-06 1.038 0.338
## Residuals 8 4.154e-05 5.193e-06
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000010 0.0000101 0.017 0.8992
## df[[3]] 1 0.003109 0.0031092 5.329 0.0543 .
## Residuals 7 0.004084 0.0005834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.0302126 1 51.7844 0.0001781 ***
## df[[1]] 0.0009266 1 1.5883 0.2479580
## df[[3]] 0.0031092 1 5.3292 0.0543059 .
## Residuals 0.0040840 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
dark.RQ <- dark %>% dplyr::select(contains("RQ"))
dark.RQ.SPF <- dark.RQ %>% dplyr::select(contains("SPF"))
dark.RQ.GF <- dark.RQ %>% dplyr::select(contains("GF"))
# Only taking last 5 KO mice
for (name in colnames(dark.RQ.SPF)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
dark.RQ.SPF %<>% dplyr::select(-name)
}
}
rcf.SPF <- lapply(dark.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(dark.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.905230988953924"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.884114999600456"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.3281, num df = 4, denom df = 4, p-value = 0.79
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1382837 12.7562514
## sample estimates:
## ratio of variances
## 1.328149
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = 2.2069, df = 7.8442, p-value = 0.05902
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.001406623 0.059374039
## sample estimates:
## mean of x mean of y
## 0.9079685 0.8789848
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2052 -0.4948 -0.2523 0.2112 1.6742
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.19 11.13 6.933 0.00012 ***
## df[[2]] -55.61 12.46 -4.464 0.00210 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.928 on 8 degrees of freedom
## Multiple R-squared: 0.7135, Adjusted R-squared: 0.6777
## F-statistic: 19.93 on 1 and 8 DF, p-value: 0.0021
##
## y = 77 + -56x, r^2 = 0.714
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.0000000000 0.0000000000 0.0000000000 0.8572243491 0.8966261382 0.0394017891
## sum median mean SE.mean CI.mean.0.95 var
## 4.3949240720 0.8868342514 0.8789848144 0.0086072838 0.0238976511 0.0003704267
## std.dev coef.var
## 0.0192464717 0.0218962506
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.000000000 0.000000000 0.000000000 0.878609328 0.940941881 0.062332554
## sum median mean SE.mean CI.mean.0.95 var
## 4.539842610 0.907323949 0.907968522 0.009919496 0.027540936 0.000491982
## std.dev coef.var
## 0.022180667 0.024428894
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 27.10000000 31.20000000 4.10000000
## sum median mean SE.mean CI.mean.0.95 var
## 141.40000000 27.50000000 28.28000000 0.75193085 2.08769472 2.82700000
## std.dev coef.var
## 1.68136849 0.05945433
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 5.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 133.70000000 26.90000000 26.74000000 0.58360946 1.62035962 1.70300000
## std.dev coef.var
## 1.30499042 0.04880293
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0244 0.8798
## 8
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.00210 0.0021001 4.87 0.0584 .
## Residuals 8 0.00345 0.0004312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.929 5.929 2.618 0.144
## Residuals 8 18.120 2.265
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 2.047e-05 2.047e-05 3.581 0.0951 .
## Residuals 8 4.574e-05 5.717e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.002100 0.0021001 11.24 0.0122 *
## df[[3]] 1 0.002142 0.0021421 11.47 0.0117 *
## Residuals 7 0.001308 0.0001868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.031750 1 169.9716 3.638e-06 ***
## df[[1]] 0.000282 1 1.5107 0.25874
## df[[3]] 0.002142 1 11.4673 0.01166 *
## Residuals 0.001308 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.0116573479895595"
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] WRS2_1.0-0 pastecs_1.3.21 multcomp_1.4-12 TH.data_1.0-10
## [5] MASS_7.3-51.5 survival_3.1-11 mvtnorm_1.1-0 effects_4.1-4
## [9] compute.es_0.2-4 car_3.0-7 carData_3.0-3 nplr_0.1-7
## [13] purrr_0.3.3 scales_1.1.0 tibble_2.1.3 reshape2_1.4.3
## [17] ggpubr_0.2.5 ggplot2_3.3.0 dplyr_0.8.5 magrittr_1.5
## [21] docstring_1.0.0
##
## loaded via a namespace (and not attached):
## [1] splines_3.6.3 assertthat_0.2.1 cellranger_1.1.0 yaml_2.2.1
## [5] pillar_1.4.3 lattice_0.20-40 glue_1.3.2 digest_0.6.25
## [9] ggsignif_0.6.0 minqa_1.2.4 colorspace_1.4-1 sandwich_2.5-1
## [13] cowplot_1.0.0 htmltools_0.4.0 Matrix_1.2-18 survey_3.37
## [17] plyr_1.8.6 pkgconfig_2.0.3 haven_2.2.0 openxlsx_4.1.4
## [21] rio_0.5.16 lme4_1.1-21 mgcv_1.8-31 farver_2.0.3
## [25] withr_2.1.2 nnet_7.3-13 klippy_0.0.0.9500 cli_2.0.2
## [29] crayon_1.3.4 readxl_1.3.1 evaluate_0.14 fansi_0.4.1
## [33] nlme_3.1-144 forcats_0.5.0 xml2_1.2.5 foreign_0.8-75
## [37] tools_3.6.3 data.table_1.12.8 hms_0.5.3 mitools_2.4
## [41] lifecycle_0.2.0 stringr_1.4.0 munsell_0.5.0 zip_2.0.4
## [45] compiler_3.6.3 rlang_0.4.5 grid_3.6.3 nloptr_1.2.2.1
## [49] labeling_0.3 rmarkdown_2.1 boot_1.3-24 gtable_0.3.0
## [53] codetools_0.2-16 abind_1.4-5 DBI_1.1.0 reshape_0.8.8
## [57] roxygen2_7.1.0 curl_4.3 R6_2.4.1 zoo_1.8-7
## [61] knitr_1.28 mc2d_0.1-18 stringi_1.4.6 Rcpp_1.0.4.6
## [65] vctrs_0.2.4 tidyselect_1.0.0 xfun_0.12